#1. 对于RNAseq_Data数据,65岁以前作为中年组，65岁以后作为老年组，使用DESeq2进行2组的差异分析。
# 加载相关的R包
library(writexl)
BiocManager::install("tidyverse")
library(tidyverse)
library(DESeq2)
library(SummarizedExperiment)

# 确保DESeq2正确加载
if(!requireNamespace("DESeq2", quietly = TRUE)) {
  BiocManager::install("DESeq2")
}
library(DESeq2)

# 设置文件地址
setwd("C:\\Users\\32076\\Desktop\\RNAseq_Data")

# 读取文件
clin_inf <- read.table("clin_inf.csv", header = TRUE)
count <- read.table("count.csv", header = TRUE)

# 处理文件数据
colnames(count) <- sub(".$", "", colnames(count)) # 把样本名的最后一个字母删掉
colnames(count) <- sub("gene_i", "gene_id", colnames(count)) # 这步用处不大，好看点
colnames(count) <- sub("X(...$)", "PTB\\1", colnames(count)) # X123改成PTB123
colnames(count) <- sub("X", "XYA", colnames(count)) # X1231244改成XYA1231244
count <- count[, order(colnames(count))] # 统一排列顺序

# 确保计数数据是整数
count[, 2:ncol(count)] <- round(count[, 2:ncol(count)])
count[, 2:ncol(count)] <- lapply(count[, 2:ncol(count)], as.integer)

# 增加一组agegroup包含old和young
clin_inf <- mutate(clin_inf, 
                   AgeGroup = ifelse(年龄 < 65, "Young", "Old"),
                   AgeGroup = factor(AgeGroup, levels = c("Young", "Old")))

# 提取计数数据
count_data <- count %>%
  column_to_rownames("gene_id") %>%
  as.matrix()

# 验证并确保样本顺序与临床信息一致
cat("样本匹配验证：\n")
matched_samples <- match(colnames(count_data), clin_inf$样本名称)
if(any(is.na(matched_samples))) {
  cat("警告：发现不匹配的样本！\n")
  cat("计数数据中存在但临床信息中不存在的样本：\n")
  print(colnames(count_data)[is.na(matched_samples)])
} else {
  cat("所有样本匹配成功！\n")
}

# 确保临床信息顺序与计数数据一致
clin_inf <- clin_inf[matched_samples, ]

# 构建DESeq2对象
dds <- DESeqDataSetFromMatrix(
  countData = count_data,
  colData = clin_inf,
  design = ~ AgeGroup
)

# 过滤低表达基因（可选但推荐）
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]

# 运行DESeq2分析
dds <- DESeq(dds)

# 提取结果（Old 和 Young）
res <- results(dds, contrast = c("AgeGroup", "Old", "Young")) %>%
  as.data.frame() %>%
  rownames_to_column("gene_id") %>%
  arrange(padj)

# 保存结果
write.csv(res, "DESeq2_AgeGroup_Results.csv", row.names = FALSE)
